The Doppler Flip in HD 100546 as a Disk Eruption: The Elephant in the Room of Kinematic Protoplanet Searches

The first exoplanets were detected indirectly through the stellar reflex motion (Mayor & Queloz 1995). Likewise, the gas velocity field in protoplanetary disks bears the imprint of planet/disk interactions, and molecular-line observations have emerged as an alternative approach to the identification of protoplanets. The first observability predictions, by Perez et al. (2015), linked embedded planets to kinks or wiggle-shaped deviations from sub-Keplerian rotation in channel maps (sky images in position–velocity data cubes, hereafter I_v ( x )). These predictions appear to match the data in HD 163296 and in HD 97048, where the most conspicuous deviations are interpreted in terms of the location and mass of the perturber through comparison with hydro simulations (Pinte et al. 2018a, 2019).

New hopes for protoplanet detection are thus coming from disk kinematics, with the interpretation of molecular-line data using hydrodynamical simulations. However, the simulations required to link the channel-map data to embedded bodies tend to idealize the systems. The unicity of a given planetary system configuration remains to be explored. It turns out that wiggles or kinks are almost ubiquitous within a disk, and also very frequently observed, as they are picked up in at least nine disks (Pinte et al. 2020) from the DSHARP ALMA Large Program (Andrews et al. 2018), despite these observations not being designed to trace small velocity structures. It may be that the ubiquity of velocity deviations reflects the global impact of gravitational interactions with the embedded planetary systems or even large-scale perturbations in gravitationally unstable disks (Hall et al. 2020; Paneque-Carreño et al. 2021). Another uncertainty is that a fraction of these structures stems from synthesis imaging limitations or from structure in the underlying optical depth. In any case, although velocity structures are expected from planet/disk interactions, it is not possible to unequivocally infer the location of embedded planets directly from the observed kinks or wiggles.

An alternative to pick up and locate embedded protoplanets has been proposed based on the velocity centroid maps, v_◦( x ) = 〈v I_v 〉 = ∫dv v I_v /∫I_v dv. After the subtraction of the axially symmetric background flow, the velocity reversal along the planetary wakes, right at the location of the planet, should be observable as a Doppler flip in molecular-line maps (Pérez et al. 2018; Casassus & Pérez 2019), i.e., as an abrupt and localized change of sign in the velocity centroid of molecular-line maps.

Very deep and long-baseline ALMA observations of HD 100546 (up to 12 km; Pérez et al. 2020) revealed the largest known velocity deviations seen in channel maps. After subtraction of the background axially symmetric flow, these deviations translate into a conspicuous Doppler flip, whose amplitude reaches one-third of the Keplerian velocity (Casassus & Pérez 2019). The root of the Doppler flip is constrained within ∼5 au, at a stellocentric radius R_f ∼ 28 ± 3 au. The amplitude of the flip requires a rather massive body, ≳ 5 M_jup for an interpretation in terms of an embedded planet. However, as noted by Casassus & Pérez (2019), its coincidence with the massive dust ring of HD 100546's disk is in contradiction with a massive protogiant, as theory predicts the clearing of a gap by the protoplanet. The only Doppler flip reported in the literature, corresponding to the largest velocity deviations known, is thus the elephant in the room of kinematic protoplanet searches. It is pointless to search for planets causing small velocity deviations, if the largest velocity deviation known remains unexplained.

In this Letter, we present new ALMA observations of HD 100546 in ¹²CO(2–1) and adjacent continuum, at unprecedented angular resolutions (Section 2). We analyze the line kinematics in terms of 3D velocity deviations from an axially symmetric flow (Section 3) and discuss these results in terms of the possible origins for the kinematic structures (Section 4) before summarizing our results (Section 5).

The data presented in this Letter correspond to ALMA observations of HD 100546, carried out in program 2018.1.01309.S (P.I. S. Pérez). The source was acquired in six executions blocks, for a total of 9.2 hr on source. The array was in configuration C-10, with baselines ranging from 83 to 15,238 m, and 41–45 active antennas. A log of observations is provided in Table 1. Quasars J1107–4449 and J1145–6954 were used for flux and phase calibration, respectively.

The correlator was set up to provide four spectral windows (spws). The continuum was sampled with spw 0, 1, and 2, at 233.0 17.0 and 219.0 GHz, each covering 2 GHz and thus yielding a total bandwidth of 6 GHz. The ¹²CO(2–1) line, at a rest frequency of 230.538 GHz, was sampled in spw 3 with 158.74 m s⁻¹ channels, which we resampled to the local standard of rest into 160 m s⁻¹ channels. The data were calibrated by staff from the North America ALMA Regional Center.

In this report, we focus on the disk kinematics and continuum emission of the central regions of the HD 100546 disk, along its inner ring, which is only ∼ 02 in radius. We report exclusively on the new long-baseline data and do not concatenate with the more compact configuration previously presented in Pérez et al. (2020). Keplerian rotation spans 35 mas in the time interval of ∼2 yr between the acquisitions of the present data and that of the intermediate-length configuration. This arc is almost 3× the finest angular scales reported here (about 012 or even 009, as measured in the central source assuming that it is unresolved, see Figure 1, or as estimated with one-third of the natural-weights beam, see Cárcamo et al. 2018) and almost twice those from Pérez et al. (2020). The concatenated continuum image of the ring is thus smoothed to angular scales two to three times coarser. However, the concatenated ¹²CO(2–1) data is much more sensitive to extended emission. Its presentation is postponed to a forthcoming article on the large-scale disk kinematics.

Zoom In Zoom Out Reset image size

The data were self-calibrated using the standard procedure. Each execution block (EB) was self-calibrated independently. We set up CASA task gaincal to average whole scans (option "solint" set to "inf"). The self-calibration process provided an improvement of a factor of 1.5–2 per EB in tclean images using Briggs weights with robustness parameter r = 0.0, with typical dynamic ranges (i.e., signal-to-noise ratio) of ∼175 to ∼200.

Following the self-calibration procedure, we performed the final synthesis imaging with the uvmem package (Casassus et al. 2006; Cárcamo et al. 2018), in the same way as in Casassus et al. (2021). In the present application, uvmem produces a positive-definite model image I^m and corresponding model visibilities ${V}_{k}^{m}$ that fit the visibility data ${V}_{k}^{\circ }$ in a least-squares sense, i.e., by minimizing

$$\begin{eqnarray}&&{\chi }^{2}=\sum _{k=1}^{{N}_{\mathrm{vis}}}{\omega }_{k}| {V}_{k}^{\circ }-{V}_{k}^{m}{| }^{2},\end{eqnarray} \tag{ 1 }$$

where ω_k correspond to the visibility weights and where the sum runs over all visibility data (i.e., without gridding).

For the restoration of the continuum image, shown in Figure 1, we use Briggs weights with robustness parameter r = 0, resulting in a clean beam 0024 × 0016/31.1, where we give the beam major axis (bmaj), minor axis (bmin), and direction (bpa) in the format bmaj×bmin/bpa. This restored image can be compared to that obtained by standard packages such as tclean in casa.

The channel maps of the ¹²CO(2–1) line, shown in Appendix B, were restored with Briggs r = 1 for a beam 0038 × 0024/0.0. In the case of emission that is much more extended than the beam, the level of aliasing is reduced with masks that isolate the signal. Here we used the same strategy as in Casassus et al. (2021), with Keplerian masks obtained with the tool ¹² developed by Teague (2020).

Moments maps were extracted with double-Gaussian fits to each line of sight, using the package GMoments ¹³ (Casassus et al. 2021). The line profile for a given line of sight is modeled as ${I}_{v}={\sum }_{j=1}^{2}{I}_{j}^{A}\exp \left(-\tfrac{{\left(v-{v}_{j}^{\circ }\right)}^{2}}{{\sigma }_{j}^{2}}\right)$. The top side, which faces the observer, is approximately traced with the brighter Gaussian, with the largest amplitude. The moment maps that result from the double-Gaussian fits are summarized in Figure 2. The following analysis of disk kinematics is based on Figure 2(c). Systemic velocity was set to 5.71 km s⁻¹, based on the kinematic analysis presented below (Section 3).

Zoom In Zoom Out Reset image size

Disk rotation can be approximated with an axially symmetric velocity field $\tilde{{\boldsymbol{v}}}({\boldsymbol{r}})$ in terms of a 3D disk rotation curve $\tilde{{\boldsymbol{v}}}(R)$, with $\tilde{{\boldsymbol{v}}}({\boldsymbol{r}})=({\tilde{v}}_{R}(R),{\tilde{v}}_{\phi }(R),{\tilde{v}}_{z}(R))$ in cylindrical coordinates. The deviations from the axially symmetric flow are then ${\boldsymbol{u}}({\boldsymbol{r}})={\boldsymbol{v}}({\boldsymbol{r}})-\tilde{{\boldsymbol{v}}}({\boldsymbol{r}})$, where v ( r ) is the 3D velocity field of the gas. If the emission stems from a thin layer at the unit-opacity surface, the first moment of the axially symmetric flow along a line of sight $\hat{s}({\boldsymbol{x}})$ is

$$\begin{eqnarray}&&{v}_{\circ }^{m}({\boldsymbol{x}})=\hat{s}({\boldsymbol{x}})\cdot \tilde{{\boldsymbol{v}}}\left({{\boldsymbol{f}}}_{\mathrm{PA},i,\psi }^{-1}({\boldsymbol{x}})\right),\end{eqnarray} \tag{ 2 }$$

where x = f _PA,i,ψ (R, ϕ). The coordinate transform f depends on the disk position angle PA and inclination i and relates x to the point, seen in the disk's cylindrical system, at the intersection of the line of sight with the disk surface, which is approximated locally as a cone with opening angle ψ. The observed velocity deviations are thus $\delta {v}_{\circ }({\boldsymbol{x}})={v}_{\circ }-{v}_{\circ }^{m}$ and derive from the intrinsic velocity deviations u ( r ).

Here we apply the ConeRot package, ¹⁴ described in Casassus & Pérez (2019) and Casassus et al. (2021), to derive the 3D rotation curve $\tilde{{\boldsymbol{v}}}(R)=({\tilde{v}}_{R}(R),{\tilde{v}}_{\phi }(R),{\tilde{v}}_{z}(R))$, in disk-centered cylindrical coordinates. Details on the extraction of this 3D rotation curve are given in Appendix C. We focus on the vicinity of the Doppler flip and postpone an analysis of the whole disk kinematics to a forthcoming article. Although ConeRot can also estimate the disk orientation based on the line data alone, these orientation estimates are partly degenerate with the rotation curve. We therefore fixed the disk orientation to that obtained from the continuum, which was extracted using the MPolarMaps package (git@github.com:simoncasassus/MPolarMaps.git; Casassus et al. 2021). This package minimizes the variance in the radial profile of the continuum intensity to estimate PA and i. The result is PA = 323 ± 07, i = 41.7 ± 04. The ring center was estimated to be offset from the phase center (itself pointed at the star) by Δα =0008 ± 00008 and Δδ = 0001 ± 00007, and the inferred systemic velocity is ${v}_{{\rm{s}}{\rm{y}}{\rm{s}}}=5.71\pm 0.02\,{\rm{k}}{\rm{m}}\,{{\rm{s}}}^{-1}$

The rotation curve for the HD 100546 disk, with meridional flows, is shown in Figure 3. The rotation curve for a purely azimuthal flow is very similar to Figure 3(a), although with a somewhat higher stellar mass (2.13 ± 0.04 < M_⋆/M_⊙ < 2.26 ±0.04). The range of possible values for the stellar mass, including both cases with and without meridional flows, is 2.12 ± 0.04 < M_⋆/M_⊙ < 2.26 ± 0.04, where the upper limits stem from vertical Keplerian shear, while the lower limit stems from cylindrical rotation (see Appendix C). These mass estimates, which assume a distance of 108.1 ± 0.4 pc (Gaia Collaboration et al. 2018), bracket the photospheric estimate of ${M}_{\star }={2.18}_{-0.17}^{+0.02}\,{M}_{\odot }$ (Wichittanakom et al. 2020). The lower-mass estimate of 1.83 ±0.01 M_⊙ reported in Casassus & Pérez (2019) was affected by the extraction of the moment map using a single Gaussian, which is more contaminated by the back side of the disk than the two Gaussians and which led to higher inclinations (of ∼ 46°; the stellar mass is sensitive to the ∼2nd power of inclination, Casassus et al. 2021, their Equation (10)), along with a missing $\sin (i)$ correction.

Zoom In Zoom Out Reset image size

Given the disk rotation curve, we now turn to the nonaxial kinematics, i.e., to local deviations from axially symmetric flow. Figure 4 presents a summary of the kinematics, as seen on the sky, and compares them with the continuum emission. The same kinematic structures are shown in a face-on deprojection in Figure 5. The deprojection technique accounts for the thickness of the CO(2–1) surface and is described in Casassus & Pérez (2019).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The velocity deviations from a purely azimuthal flow, inferred from Figures 4(c) and 5, are very similar to those reported in Casassus & Pérez (2019). There is a coherent velocity structure in the form of an arc that appears to coincide with the bright submillimeter continuum ring. The correspondence with structures in the continuum is shown in more detail in Figure 7, where we also see that the center of the flip coincides with a bifurcation along this ridge.

The Doppler flip extends over ∼120° or ∼56 au in azimuth, at a fixed orbital radius. This is best seen in the polar expansion of Figure 6, which also shows that the intrinsic velocity deviations u must reach ±1/3 Keplerian. As noted by Casassus & Pérez (2019), along the disk minor axis, the line of sight of the azimuthal components cancels, ${u}_{\phi }\hat{\phi }\cdot \hat{s}=0$. This implies that u_r or u_z contribute significantly to δ v_◦( x ) in the blue side of the flip, at least along the disk minor axis. Likewise, along the disk major axis, the line-of-sight component of a radial component cancels, ${u}_{R}\hat{R}\cdot \hat{s}=0$, which implies that the red side of the flip is either due to u_z < 0 (receding toward the back side of the disk), or to u_ϕ > 0, accelerating beyond Keplerian rotation at the ascending node.

Zoom In Zoom Out Reset image size

With an extension of the rotation curve to 3D, including meridional flows, the blue side of the flip essentially disappears from the velocity deviation map of Figure 4(d), while the red side remains. The above qualitative interpretations are thus confirmed on average, not just along the disk axes.

Interestingly, in the face-on view of the nonaxial kinematics shown in Figure 5(d), the blue and red sides of the flip align closely with the brightest continuum ridge. The velocity deviations appear to follow the peak continuum emission. Given the hints for an optically thick 225 GHz continuum (see Pérez et al. 2020; along with the reduction of peak amplitude ${I}_{1}^{A}$ and the outline of the ring in absorption in ${I}_{2}^{A}$), the bright ridge is likely a temperature feature rather than a local increase in density. Thus, the strongest velocity deviations coincide with the hottest continuum emission.

In what follows we consider the origin for the observed kinematics, taking into account the constraints on the intrinsic velocity field u . We do not consider the possibility of planet–disk interactions driven by planets inside the cavity, because the largest kinematic signatures of such bodies is in their immediate vicinity. The possibility of stellar companions inside the cavity is ruled out by the sparse-aperture-masking data presented in Pérez et al. (2020).

4.1. A Protoplanet Outflow

4.2. A Disk Eruption

Perhaps the Doppler flip could be more accurately described as a disk eruption, a surface disturbance driven by an embedded protoplanetary outflow. This outflow would not puncture the ¹²CO layer and would only push molecular material both vertically and in the direction of rotation. In this hypothesis, the planetary outflow does not escape the disk, and instead deposits momentum in disk material that is predominantly in Keplerian rotation, which thus curves around in its orbit in a fashion that is reminiscent of high-altitude volcanic plumes.

We now consider what kind of outflow would be needed to drive the observed intrinsic velocity deviations through such a disk eruption. We hypothesize that the observed δ v_◦ are compatible with the flow resulting from an expansion impulse Γ but that this flow does not puncture the surface. This impulse corresponds to the force per unit volume and mass caused by the ram pressure $\rho {v}_{s}^{2}$ inside an expanding shell with mass density ρ, width Δr_s and radius r_s ,

$$\begin{eqnarray}&&{\rm{\Gamma }}=\,\displaystyle \frac{{v}_{s}^{2}}{{\rm{\Delta }}{r}_{s}}\gtrsim \,\displaystyle \frac{{v}_{s}^{2}}{{r}_{s}}.\end{eqnarray} \tag{ 3 }$$

A shell with mass M_s will not puncture the surface if its impulse is balanced by the external pressure $\rho {c}_{s}^{2}$,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left({M}_{s}{v}_{s}\right)\sim \rho {c}_{s}^{2}4\pi {r}_{s}^{2}.\end{eqnarray} \tag{ 4 }$$

We avoid the question of the launching mechanism, and consider a shell blown at constant v_s by an embedded wind or outflow, sustained with a mass-loss rate ${\dot{M}}_{s}$:

$$\begin{eqnarray}&&{\dot{M}}_{s}{v}_{s}\sim \rho {c}_{s}^{2}4\pi {r}_{s}^{2}.\end{eqnarray} \tag{ 5 }$$

If we evaluate ρ at the CO(2–1) unit-opacity surface H₁ with a hydrostatic density profile,

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{M}}_{s}}{{{\rm{\Sigma }}}_{g}}\sim \displaystyle \frac{{r}_{s}^{2}}{{v}_{s}}\displaystyle \frac{{H}_{1}}{{f}_{1}}{e}^{-\tfrac{1}{2}{f}_{1}^{2}}{{\rm{\Omega }}}_{K}^{2},\end{eqnarray} \tag{ 6 }$$

where f₁ = H₁/H ∼ 2−5, c_s = Ω_K H, and Σ_g is the gas surface density. For the outflow to have an appreciable impact, we require r_s ≳ H. At the center of the flip in the disk of HD 100546 (M_⋆ ∼ 2 M_⊙), we have a cylindrical radius R_f ∼ 28 au and a disk aspect ratio given by ConeRot of $h=\tan (\psi )=0.2$ for ψ = 11° (which gives a midplane temperature T_m = 44 K if f₁ = 4, for a mean molecular weight μ = 2.17). If we set v_s ∼ v_f ∼ 1 km s⁻¹, then $\tfrac{{\dot{M}}_{s}}{{{\rm{\Sigma }}}_{g}}\sim {10}^{14}\,$ cm² s⁻¹. In other words, for Σ_g ∼ 5 g cm⁻², an outflow with ${\dot{M}}_{s}\sim {10}^{-6}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$ should have an appreciable impact.

If we adopt a standard ejection to accretion mass ratio of ∼10 (Konigl & Pudritz 2000), we see that the growth of a super-Earth over 1 Myr could drive such a disk eruption at an altitude over the disk midplane of 4H. It is interesting to note that the above argument, if extended to f₁ = H₁/H = 5 (so with a midplane temperature of 28 K), yields a much smaller mass-loss rate of ${\dot{M}}_{s}\sim {10}^{-8}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$, which would correspond to the accretion of a few percent of a super-Earth mass over 1 Myr. We conclude from this argument that while embedded outflows would have little impact in the denser midplane, they should lead to conspicuous velocity deviations at higher altitudes, which is consistent with the basic properties of the velocity deviations seen in the disk of HD 100546.

4.3. Disk Mass-Loading Hot Spot?

In this letter we presented new observations of HD 100546 in ¹²CO(2–1), and adjacent continuum, at unprecedented angular resolutions, which we interpreted in terms of 3D velocity deviations around the axially symmetric flow driven by Keplerian rotation. Our main conclusions are:

• 1.  
  The previously reported Doppler flip stands out as a large velocity deviation from a purely tangential and axially symmetric flow, of ±1.9 km s⁻¹, corresponding to intrinsic velocity deviations of ±1/3 Keplerian (Figure 6)
• 2.  
  However, the blue side of the Doppler flip disappears when accounting for vertical and radial flows in the axially symmetric background flow (Figures 4 or 5). Given the length of the blue side, which extends 90° in azimuth along a fairly constant radius, the bulk of the velocity deviations should thus be vertical (because the height of the emitting surface is fixed at the τ = 1 surface).
• 3.  
  The nonaxial velocity deviations are very closely aligned with the brightest ridge along the continuum ring, and the center of the flip coincides with a bifurcation along this ridge (see Figure 7).
• 4.  
  The central mass is 2.1 < M_⋆/M_⊙ < 2.3, in agreement with photospheric mass estimates (Figure 3).

Zoom In Zoom Out Reset image size

Pending dedicated hydrodynamical simulations, we tentatively propose a qualitative interpretation for the Doppler flip in terms of a surface disturbance to the Keplerian flow, i.e., a disk eruption, driven by an embedded outflow launched by a ∼10 M_Earth body. Such an outflow would deposit vertical momentum in disk material that is predominantly in Keplerian rotation and would pollute the disk material in sulfur-bearing species downstream from the source of the outflow, thus accounting for the SO line data. Another possible interpretation could involve a disk accretion hot spot at the convergence of a disk-mass-loading streamer infalling from the envelope.

Whichever the nature of the strong velocity deviations in the disk of HD 100546, the present data and analysis show that protoplanet searches based on disk kinematics should consider alternative interpretations to the non-Keplerian structures observed in molecular-line channel maps. The largest velocity deviations can be vertical or radial and seen on the disk surface, rather than azimuthal and in the disk midplane, as would be expected from the gravitational influence of an embedded and massive planet on its immediate vicinity.

We thank the anonymous referees for comments that improved the manuscript, and Richard Teague for interesting comments and suggesting Figure 10 (in the Appendix) S.C., M.C., and P.W. acknowledge support from Agencia Nacional de Investigación y Desarrollo de Chile (ANID) given by FONDECYT Regular grants 1211496, ANID PFCHA/DOCTORADO BECAS CHILE/2018-72190574, ANID project Data Observatory Foundation DO210001, ALMA-ANID postdoctoral fellowship 31180050, and FONDECYT Postdoctorado grant 3220399. This work was funded by ANID—Millennium Science Initiative Program—Center Code NCN2021_080. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01309.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Facility: ALMA. -

Software: MPolarMaps (https://github.com/simoncasassus/MPolarMaps, Casassus et al. 2021), uvmem (https://github.com/miguelcarcamov/gpuvmem, Cárcamo et al. 2018), ConeRot (https://github.com/simoncasassus/ConeRot, updated version described in this article).

The observed channel maps for the ¹²CO(2–1) data are shown in Figure 8.

Zoom In Zoom Out Reset image size

The extraction of the disk rotation curve from the velocity centroid map rests on the assumption that it corresponds to the front CO layer, as we need to convert each line of sight to a unique location in disk-centric coordinates. As explained in Section 2, we approximately separate the two τ = 1 CO(2–1) layers with a double Gaussian, where the brighter Gaussian traces the front side. While this approximation breaks down in some lines of sight, such as where the back side of the disk is exposed on the near side of the disk minor axis, a double-Gaussian fit is in any case a better representation of the observed profile than that obtained with a single Gaussian. In addition, in lines of sights where the separation of the front and back CO layers is larger than the line width (see Figure 3 of Pinte et al. 2018b, for an illustration of the two CO layers), the line profile is double peaked, which for a single-Gaussian fit would bias the velocity centroid in between the two velocity components.

The separation of the two CO(2–1) layers results in qualitatively different structures for ${v}_{1}^{\circ }$ and ${v}_{2}^{\circ }$. The position of the peak velocity centroid of the brighter Gaussian (${v}_{1}^{\circ }$), in absolute value, drifts away from the disk major axis with distance from the star, roughly toward the northeast, as expected for a flared surface (see Figure 2(c)). The opposite holds for the fainter Gaussian (${v}_{2}^{\circ }$), whose peak velocity drifts toward the southwest (see Figure 2(f)). In addition, the continuum ring appears to coincide with a region of lower secondary Gaussian amplitude (${I}_{2}^{A}$; see Figure 2(e)) and can be outlined (roughly) in absorption. The separation of the two CO(2–1) layers, using the double-Gaussian fits, is illustrated in selected channels in Figure 9.

Example line of sights spectra are shown in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The continuum ring in the disk of HD 100546 has a rather small radius, of about 5 to 10 beams even in these long-baseline data, and the central cavity is devoid of CO(2–1) emission. In addition, the standalone long-baseline data of HD 100546 are not ideal to trace the larger scales. As a proof of concept for the use of the double-Gaussian fit to separate the two Gaussian layers, we have applied GMoments to archival data of HD 163296 in ¹²CO(2–1) (from the DSHARP ALMA Large Program (Isella et al. 2018), which we reimaged in the same way as for the HD 100546 data presented here). The intensity maps for each CO layer are shown for selected channels in Figure 11. The corresponding moment maps are shown in Figure 12, where the concentric continuum ring system of HD 163296 is seen in absorption in ${I}_{2}^{A}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use the same notation and reference frames as in Casassus & Pérez (2019). The sky frame is represented by ${ \mathcal S }$ and is oriented using Cartesian coordinates (x, y, z), with $\hat{y}$ aligned due north. A rotation about the vertical axis $\hat{z}$ into frame ${{ \mathcal S }}^{{\prime} }$ aligns $\hat{y}^{\prime}$ with the disk major axis, while the $(\hat{x}^{\prime} ,\hat{y}^{\prime} )$ plane is still parallel to the sky. A second rotation about axis ${\hat{y}}^{{\prime} }$, by an inclination angle i, matches the $(\hat{x}^{\prime} ,\hat{y}^{\prime} )$ plane with the disk midplane in frame ${{ \mathcal S }}^{{\prime\prime} }$.

If all of the emission originates from the top side of the disk, which is facing the observer, we have a bijection between the line of sight x and the polar coordinate of its intersection with the surface of the cone representing the disk surface, defined as the surface of unit opacity, with height H₁(R) above the midplane. This cone has opening angle $\psi =\arctan ({h}_{1})$, where h₁ = H₁/R is the aspect ratio of the disk surface. We thus transform the sky coordinates x to cylindrical polar coordinates on the surface of a cone (R, ϕ), as measured in ${{ \mathcal S }}^{{\prime\prime} }$. The coordinate transform from Cartesian coordinates in ${{ \mathcal S }}^{{\prime\prime} }$ into cylindrical coordinates in ${{ \mathcal S }}^{{\prime\prime} }$,

$$\begin{eqnarray}&&{\boldsymbol{x}}={{\boldsymbol{f}}}_{\mathrm{PA},i,\psi }(R,\phi ),\end{eqnarray} \tag{ C1 }$$

depends on disk orientation and is invertible. Similar coordinate transforms have been used in, for example, Rosenfeld et al. (2013) or in Isella et al. (2018). We implement f with the following formulae, already provided in Casassus & Pérez (2019) but which we reproduce here for completeness:

$$\begin{eqnarray}\begin{array}{rcl}{x}^{{\prime} } & = & R\sin (\phi )/\cos (i)\\ & & +({H}_{1}-R\sin (\phi )\tan (i))\sin (i),\end{array}\end{eqnarray} \tag{ C2 }$$

$$\begin{eqnarray}&&{y}^{{\prime} }=R\cos (\phi ).\end{eqnarray} \tag{ C3 }$$

The inversion of Equations (C2) and (C3) to obtain $(r,\phi )={{\boldsymbol{f}}}_{\mathrm{PA},i,+\psi }^{-1}({x}^{{\prime} },{y}^{{\prime} })$ can be achieved by noting that R is the root of

$$\begin{eqnarray}&&\displaystyle \frac{{y}^{{\prime} 2}}{{R}^{2}}+\displaystyle \frac{{\left({x}^{{\prime} }-{H}_{1}(R)\sin (i)\right)}^{2}}{{R}^{2}{\cos }^{2}(i)}=1,\end{eqnarray} \tag{ C4 }$$

and with

$$\begin{eqnarray}&&\cos (\phi )={y}^{{\prime} }/R.\end{eqnarray} \tag{ C5 }$$

If all of the observed signal stems from the top side of the disk and in the H₁(R) surface, then the first moment of the axially symmetric flow along line of sight $\hat{s}({\boldsymbol{x}})$ is, using Equation (C1),

$$\begin{eqnarray}&&{v}_{\circ }^{m}({\boldsymbol{x}})=\hat{s}({\boldsymbol{x}})\cdot \tilde{{\boldsymbol{v}}}\left({{\boldsymbol{f}}}_{\mathrm{PA},i,\psi }^{-1}({\boldsymbol{x}})\right)\end{eqnarray} \tag{ C6 }$$

or

$$\begin{eqnarray}&&\begin{array}{l}{v}_{\circ }^{m}(R,\phi ,z={H}_{1})={\tilde{v}}_{\phi }(R)\cos (\phi )\sin (i)\\ \quad +\,{\tilde{v}}_{r}(R)\sin (\phi )\sin (i)+{\tilde{v}}_{z}(R)\cos (i),\end{array}\end{eqnarray} \tag{ C7 }$$

where all coordinates are measured in ${{ \mathcal S }}^{{\prime\prime} }$, i.e., in the frame of the disk. Similar formulae for the model velocity centroid have been proposed by Teague et al. (2019).

The origin of ϕ coincides with the PA of the disk major axis, defined as the line of ascending nodes (so the red side). The sign convention used here results in ${\tilde{v}}_{r}(R)\lt 0$ for accretion. In case of prograde rotation, with i < 90°, a wind flowing away from the midplane corresponds to ${\tilde{v}}_{z}(R)\gt 0$. But for retrograde rotation, the top side of the disk that is facing the observer corresponds to z < 0 in ${{ \mathcal S }}^{{\prime\prime} }$, so the disk opening angle ψ < 0, and for a wind ${\tilde{v}}_{z}(R)\lt 0$. The observed velocity deviations are thus $\delta {v}_{\circ }({\boldsymbol{x}})={v}_{\circ }-{v}_{\circ }^{m}$ and derive from the intrinsic velocity deviations u ( r ):

$$\begin{eqnarray}\begin{array}{rcl}\delta {v}_{\circ }({\boldsymbol{x}}) & = & {u}_{\phi }(R,\phi ){| }_{z={H}_{1}}\cos (\phi )\sin (i)\\ & & +{u}_{R}(R,\phi ){| }_{z={H}_{1}}\sin (\phi )\sin (i)\\ & & +{u}_{z}(R,\phi ){| }_{z={H}_{1}}\cos (i).\end{array}\end{eqnarray} \tag{ C8 }$$

Given the disk orientation (PA, i, ψ), we solve for the disk rotation curve ${\tilde{{\boldsymbol{v}}}}_{\circ }(R)$ by performing a least-squares fit of ${v}_{\circ }^{m}(R,\phi )$ in Equation (C7) to v_◦( f _PA,i,ψ (R, ϕ)), which is the observed velocity centroid resampled by f _PA,i,ψ . We obtain ${\tilde{v}}_{\phi }({R}_{l})$, ${\tilde{v}}_{r}({R}_{l})$, and ${\tilde{v}}_{z}({R}_{l})$ for each discretized value of R_l by minimizing

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\tilde{v}}^{2} & = & \displaystyle \sum _{k\in {{ \mathcal D }}_{\phi }}w({R}_{l},{\phi }_{k})\left({v}_{\circ }({R}_{l},{\phi }_{k})\right.\\ & & {\left.-{v}_{\circ }^{m}({R}_{l},{\phi }_{k})-{v}_{s}\right)}^{2},\end{array}\end{eqnarray} \tag{ C9 }$$

which assumes that the systemic velocity v_s is known. In practice, we solve for the root of ∇ χ² = 0 using the Cramer rule. The sum in Equation (C9) runs over all azimuths within a domain ${{ \mathcal D }}_{\phi }$, with weights

$$\begin{eqnarray}&&w({R}_{l},{\phi }_{k})=\displaystyle \frac{1}{{N}_{b}({R}_{l}){\sigma }_{\circ }^{2}({R}_{l},{\phi }_{k})},\end{eqnarray} \tag{ C10 }$$

where N_b is the number of beam major axes along a circle with radius R_l , clipped so that N_b > 1. N_b approximately corrects for correlated pixels in. The domain in azimuth ${{ \mathcal D }}_{\phi }$ is typically [0, 2π], but it is interesting to restrict ${{ \mathcal D }}_{\phi }$ to [0, π], corresponding to the far side of the disk. The model velocity centroid ${v}_{\circ }^{m}(R,\phi )$ represents the axially symmetric flow and can be converted back to ${{ \mathcal S }}^{{\prime\prime} }$ by resampling with ${{\boldsymbol{f}}}_{\mathrm{PA},i,\psi }^{-1}$.

If v_s is not known, then we initially restrict to a purely azimuthal rotation curve and fit for v_s (R_l ) and ${\tilde{v}}_{\phi }({R}_{l})$ simultaneously, for all {R_l } in a relatively narrow radial domain [R₁, R₂], which we call a region. This optimization also includes the disk orientation (see next paragraph), but for constant orientation parameters (i_◦, PA_◦, ψ_◦) representative of this region. The value of v_s is then fixed to its median, and its uncertainty to its standard deviation, both extracted over overlapping regions that cover a wider range in radii. The full rotation curve is then calculated using this global v_s value.

Disks are flared, i.e., the disk aspect ratio ${h}_{1}(R)\,=\tan (\psi (R))/R$ is a function of R. With molecular-line tracers, we expect that h₁(R) will initially increase with R and then drop to zero as the tracers becomes optically thin at the outer edge of the disk. Disks may also warp, so that disk orientation will depend on R. We calculate the disk orientation profile (i(R), PA(R), ψ(R)) by performing the optimization in overlapping radial bins. The full radial extension of the disk is divided into M overlapping radial bins ${\{[{R}_{1j},{R}_{2j}]\}}_{j=1}^{M}$, thus defining radial regions $\{{{\rm{\Theta }}}_{j}({{\boldsymbol{x}}}^{{\prime} })\}{}_{j\,=\,1}^{M}$, where ${{\rm{\Theta }}}_{j}({{\boldsymbol{x}}}^{{\prime} })=1$ if $R={f}_{r}^{-1}({{\boldsymbol{x}}}^{{\prime} })\in [{r}_{1j},{r}_{2j}]$, and ${{\rm{\Theta }}}_{j}({{\boldsymbol{x}}}^{{\prime} })=0$ otherwise. In each radial bin, we maximize the log-likelihood function, −0.5χ² to obtain (i_j , PA_j , ψ_j ) with

$$\begin{eqnarray}&&\begin{array}{l}{\chi }_{\mathrm{var}}^{2}=\displaystyle \frac{1}{{\sigma }_{\mathrm{ref}}^{2}}\displaystyle \sum _{l={l}_{1}}^{{l}_{2}}\\ \quad \times \displaystyle \frac{{\sum }_{k=0}^{{N}_{\phi }-1}w({R}_{l},{\phi }_{k}){\left({v}_{\circ }({R}_{l},{\phi }_{k})-{v}_{\circ }^{m}({R}_{l},{\phi }_{k})\right)}^{2}}{{\sum }_{k=0}^{{N}_{\phi }-1}w({R}_{l},{\phi }_{k})},\end{array}\end{eqnarray} \tag{ C11 }$$

which is the radial sum of the weighted azimuthal variance of residuals relative to the expected variance ${\sigma }_{\mathrm{ref}}^{2}$, with σ_ref typically of order 0.1 km s⁻¹ when including image synthesis and channelization systematics. The sums extend over an interval in radius [R₁, R₂] and over all azimuths. We thus obtain an orientation profile ${\left\{i({R}_{j}),\mathrm{PA}({R}_{j}),\psi ({R}_{j})\right\}}_{j=1}^{M}$ as well as axially symmetric models $\{{v}_{j}^{m}({\boldsymbol{x}})\}{}_{j=1}^{M}$ in all radial regions. The axially symmetric model corresponding to these profiles is approximated by the average over all regions, i.e.,

$$\begin{eqnarray}&&{v}_{\circ R}^{m}({\boldsymbol{x}})=\displaystyle \frac{{\sum }_{j=1}^{M}{v}_{\circ }^{m}({\boldsymbol{x}}){{\rm{\Theta }}}_{j}({\boldsymbol{x}})}{{\sum }_{j=1}^{M}{{\rm{\Theta }}}_{j}({\boldsymbol{x}})}.\end{eqnarray} \tag{ C12 }$$

This approximation ignores shadowing between regions along a line of sight, and is thus only applicable at low to moderate inclinations.

The minimization of Equation (C11) and the optimal orientation parameters, along with their associated errors, were achieved with a Markov Chain Monte Carlo ensemble sampler (Goodman & Weare 2010). We use the emcee package (Foreman-Mackey et al. 2013), with flat priors and with typically 300 iterations and 30 walkers.

The azimuthal rotation curve ${\tilde{v}}_{\phi }$ can be used to constrain the stellar mass under the assumption of pure Keplerian rotation. As noted in Casassus & Pérez (2019), the stellar mass is then bracketed by two extreme cases. Assuming a rigid vertical structure, in which the midplane velocity is the same as that at height H₁, sets a lower limit. Pure vertical Keplerian shear, ignoring radial hydrostatic support, sets an upper limit, in which midplane azimuthal velocities are ${\tilde{v}}_{\phi }(R){\left(1+{h}_{1}^{2}\right)}^{3/4}$.